Electrophoretic Properties of Highly Charged Colloids: A Hybrid MD/LB Simulation 

Study 
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Using computer simulations, the electrophoretic motion of a positively charged colloid (macroion) 
in an electrolyte solution is studied in the framework of the primitive model. In this model, the 
electrolyte is considered as a system of negatively and positively charged microions (counterions 
and coions, respectively) that are immersed into a structureless medium. Hydrodynamic interac- 
tions are fully taken into account by applying a hybrid simulation scheme, where the charged ions 
(i.e. macroion and electrolyte), propagated via molecular dynamics (MD), are coupled to a Lattice 
Boltzmann (LB) fluid. In a recent electrophoretic experiment by Martin-Molina et al. [J. Phys. 
Chem. B 106, 6881 (2002)], it was shown that, for multivalent salt ions, the mobility n initially 
increases with charge density a, reaches a maximum and then decreases with further increase of a. 
The aim of the present work is to elucidate the behaviour of /i at high values of a. Even for the 
case of monovalent microions, we find a decrease of fi with a. A dynamic Stern layer is defined 
that includes all the counterions that move with the macroion while subject to an external electrical 
field. The number of counterions in the Stern layer, go, is a crucial parameter for the behavior of {i 
at high values of a. In this case, the mobility fi depends primarily on the ratio qo/Q (with Q the 
valency of the macroion). The previous contention that the increase in the distortion of the electric 
double layer (EDL) with increasing a leads to the lowering of \x does not hold for high a. In fact, 
we show that the deformation of the EDL decreases with increase of a. The role of hydrodynamic 
interactions is inferred from direct comparisons to Langevin simulations where the coupling to the 
LB fluid is switched off. Moreover, systems with divalent counterions are considered. In this case, 
at high values of a the phenomenon of charge inversion is found. 

PACS numbers: 75.60.Ej,75.50.Lk 
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I. INTRODUCTION 

The electrophoretic motion of charged macroions in an 
electric field E is a longstanding issue in colloid science 
0, 1 I, S, !, I, 0, 1 1 M, EH ■ A complex interplay oc- 
curs between hydrodynamic, electrostatic, and thermal 
forces. A macroion of charge Q drags along an electric 
double layer (EDL) of small counterions with it and the 
number of charges "bound" in the EDL depends on the 
charge density, surface properties, ion concentration and 
liphophility of the macroion as well as the specific prop- 
erties of counterions and salt ions [l2j . Since it is difficult 
to experimentally probe and control these different pa- 
rameters at microscopic length and time scales, the inter- 
pretation of experimental results often remains unclear. 
There is, e.g., no clear understanding of the variation of 
the mobility // with the surface charge density a of the 
macroion [llElEl- 

Most of the theoretical studies on electrokinetic phe- 
nomena consider the Coulomb interactions between 
macroions on the level of the linearized Poisson- 
Boltzmann equation [1, 0, 0]. This equation, valid for 
weakly charged macroions and low salt concentration, de- 
scribes the formation of an EDL. It results in a screened 
Coulomb potential (Yukawa potential) as the effective in- 
teraction potential between the macroions. The central 
parameter of this potential is the Debye screening length 
Ad = 1/k (with k the screening parameter) that mea- 
sures the spatial extent of the EDL. When a charged 



colloid moves in the presence of an electric field, the 
spherical EDL is distorted and it is coupled to the hy- 
drodynamic flow of the solvent. In theoretical approaches 
d, [H, EE E3, EH , hydrodynamic effects are incorporated 
by coupling the Poisson equation for the electrostatics 
to the Navier-Stokes equations. In the resulting elec- 
trokinetic equations Q , so-called £ potential is a central 
quantity. It is defined as the electrostatic potential ip(z s ) 
at a distance z s from the colloidal surface where the fluid 
around the colloid is at rest with respect to the colloid 
motion. The imaginary surface at distance z s from the 
colloidal surface is called surface of shear or slipping sur- 
face [H [!, [I?], EH- Some of the counterions that are 
between the surface of the macroion and the surface of 
shear are often assumed to move along with the macroion 
and define the so-called dynamic Stern layer 0, E3| ■ 

The Q potential can be related to the electrophoretic 
mobility fi = VyyjE (with Vm the steady-state macroion 
velocity). In the Hclmholtz- Smoluchowski limit, i.e. for 
kRm oo (with Rm the radius of the macroion), the 
mobility of the macroion can be obtained from the £ po- 
tential via /i = eC,{z s )/rj [1, E1,EB|, where e is the permit- 
tivity of the surrounding medium and r\ is the shear vis- 
cosity of the fluid. In the opposite Huckel-Onsager limit, 
Ad >> i?M, the mobility is given by /j, = 2e((z s )/3r) 
The definition of the potential ((z s ) at the sur- 
face of shear is usually based on no-slip boundary con- 
ditions at the colloidal surface [l8|. Since kRm — > oo 
implies a very small Debye length, the distance of the 
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surface of shear from the colloidal surface would be typ- 
ically in the nanometer range. At such length scales the 
validity of no-slip boundary condition is not clear (al- 
though a hydrodynamic description might be still valid 
(l9|). Calculations of /i which cover the entire range kRm 
from the Helmholtz-Smoluchowski limit to the Hiickel- 
Onsager limit have been done by Henry [3j. But these 
calculations, valid for low £, assumed that the counterion 
charge density is unaffected by the applied field. 

O'Brien and White 0] did extensive numerical calcula- 
tions to calculate fj, versus £ for different values of kRm in 
the framework of electrokinetic equations. For kRm > 3, 
they obtained a pronounced maximum when plotting fi 
as a function of £. This behaviour is attributed to the 
competition between the driving force due to E and the 
retarding forces due to the distortion of the charge cloud 
constituting the EDL. The driving force increases lin- 
early with Q (which is proportional to Q), whereas the 
retarding force is proportional to £ 2 . Thus for high val- 
ues of \x decreases with £ because then the retarding 
force dominates. Moreover, O'Brien and White found 
that for multivalent ions, the position of the maximum 
shifts to lower values of (. They argue that the distortion 
of the double layer is increased by multivalent countcri- 
ons resulting in the increase of retardation forces against 
the motion of the colloid. Qualitatively similar results 
were also obtained by the non-linear solution by Oshina 
et al. and by perturbation expansion in powers of £ 
by Booth [§] and Overbeek [Icj . There are also other 
models of more microscopic origin which look at charged 
systems, but these models focus mainly on the structure 
and dynamics of the dynamic Stern layer in planar geom- 
etry [I?], HH rather than on electrophoretic properties of 
charged spherical colloids. Moreover, some of these stud- 
ies disregard hydrodynamic interactions [ItJ or consider 
only zero-temperature properties (20| . 

Electrophoresis experiments on colloidal suspensions 
[3 E3, ED, S m HE HI allow to determine ac- 
curately the electrophoretic mobility fi as a function of 
different parameters such as kRm, surface charge den- 
sity of colloids, etc. However, it is difficult to disentangle 
hydrodynamic effects from those due to electrostatic in- 
teraction. Moreover, in order to fit experimental data to 
theoretical models, "renormalized charges" of the colloids 
have to be introduced [25j |. 

In a recent experimental study by Martin-Molina et 
al. [HI , a maximum was found in the electrophoretic mo- 
bility //, plotted as a function of the macroion's surface 
charge density a (note that there is a linear relationship 
between a and the C potential in the framework of the lin- 
earized Poisson-Boltzmann equation). The latter maxi- 
mum occurs in the case of a 2:1 (counterion-chargexoion- 
charge) electrolyte. A similar result can be also found in 
the numerical paper by O'Brien and White Q for a 2:1 
electrolyte. However, the experiments were performed 
for systems of highly charged macroion's. In this case, 
it is not clear whether the calculations by O'Brien and 
White yield a correct description, even on a qualitative 



level. In this work, we aim at addressing this issue by us- 
ing computer simulation techniques that may elucidate 
electrophoresis with highly charged particles on a micro- 
scopic level. 

Molecular dynamics (MD) simulations of colloidal sys- 
tems suffer from the large separation in length and time 
scales of solvent and colloidal particles. The problem of 
simulating millions of solvent particles to account for hy- 
drodynamic effects can, however, be circumvented by us- 
ing simulation techniques such as the Lattice Boltzmann 
(LB) [27[ method or other Navier-Stokes equation solver 
(see, e.g. [28j) to model the solvent. In this work, we em- 
ploy a hybrid LB /MD method [29|, [3(| to investigate the 
electrophoretic motion of a highly charged macroion in 
an electrolyte solution, thereby modeling macroion and 
electrolyte in the framework of the so-called primitive 
model (in a similar way as in Ref. [3l|). Our aim in par- 
ticular is to investigate those structural and dynamical 
aspects of the EDL which are out of scope of mesoscopic 
theories. 

In the following, we are interested in the case of high 
values of a, i.e. the regime beyond the aforementioned 
maximum in the /i — er curve. In particular the coun- 
terion distribution around the macroion is measured to 
provide insight into the electrical retarding forces affect- 
ing the mobility. In order to disentangle hydrodynamic 
from electrostatic effects, the hydrodynamic medium (LB 
fluid) is also switched off, thus simulating the system via 
Langevin dynamics. Systems with monovalent and di- 
valent salt ions are considered. In the latter case, the 
phenomenon of charge inversion is observed. 

The rest of the paper is organized as follows: We de- 
scribe the LB/MD method in Sec. II. In Sec. IIIA, we 
present the results for the monovalent electrolyte solu- 
tion (microions) and in Sec. IIIB, we focus on the charge 
inversion phenomenon with divalent microions. Finally, 
conclusions are given in Sec. IV. 

II. MODEL AND DETAILS OF THE 
SIMULATION 

The electrophoretic motion of a macroion in an elec- 
trolyte solution is studied in the framework of the so- 
called primitive model. We consider a system of a 
macroion of charge Q = Zy^e (mass M = 60 a.u.) 
and microions of charge Z ct = — le (counterions) and 
Z co = le (coions), each of which are of mass 4 a.u.. The 
interaction potential between a particle of type a and a 
particle of type (3 (a,/3 = M, ct,co) separated by a dis- 
tance r from each other is given by 

= a_/3f_ +j4 a/3 exp{-S a fl(r- cr a g)} (1) 

where e is the elementary charge and e the dielectric 
constant. We choose the value e = 80e (with cq the 
vacuum dielectric constant) for water at room tempera- 
ture. The parameters a a p denote the distance between 
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two ions at contact, a a/ 3 = R a + Rp, where R a is the 
radius of an ion of type a. In the following, we use val- 
ues for i?M that vary from 10 A to 20 A. For the mi- 
croions, we set i? c t = R co = 1 A. The exponential in 
Eq. {T]) is an approximation to a hard sphere interac- 
tion for two ions at contact. For the parameters A a p we 
choose ^4mm = 1-84 eV, A^ ct = ^4m,co = 0.05565 eV, 
and A ct . ct = A cttCO = ^4 COiCO = 0.0051 eV. The pa- 
rameters B a/ 3 are all set to 4.0 A -1 . The long-ranged 
Coulomb part of the potential and the forces were com- 
puted by Ewald sums in which we chose ole — 0.05 for 
the constant and a cutoff wavenumber k c — 2tt^/66/L 
in the Fourier part [32l . |33| . The linear dimension L of 
the simulation box is L = 160 A, using periodic boundary 
conditions. All simulations were done at the temperature 
T = 297K. 

Now, the crucial step is to model the solvent. To a first 
approximation, the effect of the solvent on colloidal parti- 
cles (macroions) can be described by a Langevin equation 
where one assumes that on the typical time scale of the 
colloidal particles their collisions with solvent particles 
are due to Gaussian random forces f r ^. These forces lead 
to a systematic friction force — £oVj(i) on the colloids, 
where £o is the friction coefficient and Vj(t) the velocity 
of particle i at time t. The resulting equations of motion 
are 

M-^ = F c>i -$ ) V i (*)+f r>i) (2) 

where R^ and Vj (t) denote respectively position and ve- 
locity of a colloidal particle i (i — 1,...,N with N the to- 
tal number of colloids) and F C) j is the conservative force 
acting on the particle. The Cartesian components of the 
random forces, (a = x, y, z), are uncorrelated random 
numbers with zero mean, i.e. 

</£(Ri,*)) = (3) 

(/^(Ri,*)/^(B4,f)) = A6 aP 5(R t -RW(t-t')(4) 

The amplitude A is determined by the fluctuation- 
dissipation theorem, A = 2fcsT£o- In our model, the mi- 
croions are also Brownian particles and in thermal equi- 
librium with the solvent. 

Equations j2])-((4]) still do not provide a complete de- 
scription of the dynamics in colloidal suspensions. In 
these equations the mass and momentum transport by 
the solvent is ignored that leads to the hydrodynamic in- 
teractions between the colloidal particles. However, hy- 
drodynamic interactions can be incorporated rather eas- 
ily into the Langevin description by replacing the abso- 
lute velocity V((t) of particle i in Eq. ^ by its velocity 
relative to the fluid velocity field u(Rj, t) at the position 
of the particle R^. The "hydrodynamic force" on particle 
i is then given by 

Fhyd.i = [Vi(i) - u(Ri, t)} + f 1Vi (5) 

In a simulation, the velocity field u(Rj,i) can be calcu- 
lated from any Navier-Stokes equation solver, whereby 



thermal fluctuations have to be included in the frame- 
work of fluctuating hydrodynamics. In this work, we use 
a Lattice-Boltzmann (LB) [2?], H3] scheme to compute 
u(Ri, t). Since the LB method yields the velocity field u 
on a lattice, an interpolation scheme has to be used to 
determine u at the position of the particle R^. Very re- 
cently, Ahlrichs and Diinweg [35j have proposed a hybrid 
MD/LB scheme on the basis of the frictional coupling 
force, Eq. ([5]), to simulate polymers in solution. 

Up to now, we have considered the colloidal particles 
as point particles. However, real colloids are extended 
objects with rotational degrees of freedom. The problem 
is to "implant" an extended object such as a sphere of 
radius i?u into a LB fluid. To this end, we follow Ladd's 
method [3~ij and represent the sphere (or any other ob- 
ject) by uniformly distributed boundary points on its sur- 
face, where the surface is permeable for the LB fluid. We 
have recently proposed a sphere model with 66 bound- 
ary nodes of mass M/66 with M being the total mass 
of the colloidal particle . This model for a spherical 
particle is also considered in the following. Each of the 
boundary points is coupled to the LB fluid by the force 
given by Eq. ([5]) . The total force on the particle is deter- 
mined by the sum over the forces on the boundary nodes 
(of course, according to Newton's third law, these forces 
with opposite sign are given back to the LB fluid) . From 
this total force the torque on the particle is calculated, 
and the total force and torque are then used to update 
the translational and rotational velocity of the particle, 
respectively. More details on our MD /LB scheme can be 
found elsewhere [29| . 

For the LB fluid a standard D3Q18 model was used, 
with which we solved the linearized Navier-Stokes equa- 
tions (the details of the D3Q18 model can be found in 
review articles and the book by Succi [13, HU). We con- 
sider an incompressible fluid in the creeping flow limit 
(zero Reynolds number) in our studies. The kinematic 
viscosity was set to v = 0.0238a 2 /r with a the lattice 
constant and r the elementary time unit of the LB fluid. 
The density of the LB fluid was set to p = 1.0mo/a 3 (mo: 
mass unit of the LB fluid). Thermal fluctuations were in- 
troduced via the addition of Gaussian random numbers 
to the stress tensor as proposed by Ladd [34 ]. The LB 
fluid is modelled on a cubic lattice with 40 3 lattice nodes, 
thus the lattice constant is a — 4 A (for the linear dimen- 
sion L — 160 A of the simulation box). The counterions 
and coions are considered as point particles with respect 
to the interaction with the LB fluid. The 66 boundary 
nodes of the macroion are placed on a sphere of radius 
i?H = 10 A. The effective hydrodynamic radius and the 
viscous retarding force on the macroion is thus deter- 
mined by i?H- On the other hand, the variation of Rm 
from 10 A to 20 A (see next section) allows a change of 
the macroion's surface charge a without changing the hy- 
drodynamic coupling of the particle to the LB fluid. For 
the friction constant £o in the coupling force, Eq. ([5]), a 
total value of 6.6mo/r is assigned to the macroion which 
corresponds to £o = OAmo/r for each of the 66 bound- 
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Q [e] 


N c t 


iVco 


k [A- 1 ] 


kRm 


121 


471 


350 


0.127 


2.54 


255 


555 


300 


0.133 


2.66 


351 


651 


300 


0.137 


2.74 


401 


651 


250 


0.126 


2.52 


601 


801 


200 


0.141 


2.82 


801 


1001 


200 


0.154 


3.08 



TABLE I: Charge of the macroion Q, number of counterf- 
oils Net , number of coions 7V co , and the corresponding values 
of the screening parameter k and the dimensionless quantity 
kRm that were used in the simulations for _Rm = 20 A. 



ary nodes. In a recent publication [291 ], we have shown 
that this value of £o corresponds to nearly stick bound- 
ary conditions. For the microions, a different value for 
the friction coefficient is used, denoted by £ofc in the fol- 
lowing. Below we discuss the influence of £ofc on the elec- 
trophoretic mobility (see Fig.[7|). Unless otherwise noted, 
the value £ofc = 0.025to /t was chosen. With respect to 
the LB fluid, the microions are treated as point particles, 
but we remind the reader that their Coulomb radii are 
set to i? ct = R co = 1 A. 

We have done MD /LB simulations of a single macroion 
of charge Q in an electric field E x pointing in the pos- 
itive x direction. The charge Q was varied from 121 e 
to 801 e. The number of counterions and coions (iV ct 
and AT co , respectively) used for a given value of Q are 
listed in Table [IJ Also included in this list is the Debye 
screening parameter k = [47rAs(n c t + n co ) / 'L 3 ] 1 / 2 which 
was roughly kept constant around 0.13 A -1 (the Bjerrum 
length Ab = e 2 /(47re r eofcsT') is in our case equal to about 
7 A). For comparison, we also carried out simulations 
with a "Langevin dynamics" (LD) where the coupling to 
the LB fluid was switched off, i.e. u(R, t) = 0. For the 
LB/MD and LD simulations, the equations of motion 
were integrated by a Heun algorithm with a time step of 
1 fs. This very small time step is necessary because we 
consider explicitly counterions and coions as microscopic 
particles. 

The simulations were done as follows: We first equili- 
brated the system for 25000 time steps without electric 
field and without coupling to the LB fluid. Then, the 
system was coupled to the LB fluid and the electric field 
was switched on, followed by simulation runs over 400000 
time steps. After 100000 time steps, the steady state was 
reached and the positions and velocities of the ions were 
stored every 500 time steps to determine the averaged 
quantities such as the steady state velocity Vm of the 
macroion. 



III. RESULTS 

Now the simulation results are presented for the 
steady-state electrophoretic motion of a macroion in an 
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FIG. 1: a) V c t(r)/VM as a function of r/RM for different val- 
ues of the electric field, as indicated. Data is shown for 1 
macroion of charge Q = 255e, Net = 555 monovalent coun- 
terions and iVco = 300 monovalent coions. b) Cumulative 
counterion charge q(r) for the same parameters as in a), go is 
the cumulative counterion charge at which V c t (r) is zero. The 
determination of qo is indicated in the figure by the dashed 
lines for the example E = 0.01 V/A. 



electrolyte. Systems of highly charged macroions are 
considered, i.e. their surface charge densities vary be- 
tween a = 0.02 e/A 2 (38 ^C/cm 2 ) and a = 0.3 e/A 2 
(500 ^C/cm 2 ). The density of coions (iV co = 300) in 
our simulation box corresponds to a salt concentration 
of 0.012 mol/1. 



A. Monovalent microions 

In this section, we consider electrolyte solutions that 
consist of monovalent counter- and coions. For a micro- 
scopic understanding of electrophoresis, it is of particular 
interest to display the dynamic distribution of counte- 
rions in the vicinity of the macroion, i.e. in the EDL. 
This distribution is determined by an interplay between 
the electrostatic attraction and the hydrodynamic flow 
around the macroion. Some of the counterions in the 
EDL move along the same direction as the macroion, 
whereas, due to the electric field, other counterions in 
the EDL are accelerated in the opposite direction. We 
analyze the dynamic behavior of the counterions by their 
average velocity V c t (r) as a function of the radial distance 
r from the center of the macroion and the cumulative 
counterion charge q(r) around the macroion. 

First, we consider a system consisting of one macroion 
of charge Q = 255e, N ct — 555 monovalent counterions 
and N co = 300 monovalent coions. For this case, Fig. [T] 
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displays V c t{r)/V~M (with Vm the steady state macroion 
velocity) and q{r) for different choices of the electric 
field E. Within the statistical errors, the results essen- 
tially coincide for the two lowest values of the electric 
field, E = 0.01 V/A and E = 0.008 V/A. The choice 
E = 0.01V/A is used for all further results presented 
in this work, in particular for the estimates of the elec- 
trophoretic mobility, fi = Vm/E. As shown previously 
[29j |. for Q > 255e the linear response regime is essen- 
tially achieved for E < 0.01 V/A. 

From the behavior of V c t(r)/VM, different regions can 
be identified with respect to the distance from the cen- 
ter of the macroion. Close to the macroion's surface a 
layer with a thickness of about 0.25 Rm is formed where 
Vct/VM is constant with a value around 0.8 for the two 
lowest values of E. In the following, this region is called 
the dynamic Stern layer, where, due to the electrostatic 
attraction by the macroion, counterions are essentially 
condensed onto the surface of the macroion. Beyond the 
Stern layer, the counterion velocity changes quickly its 
sign, thus indicating a motion in the direction opposite 
to that of the macroion. The point where the counterion 
velocity vanishes can be used as a measure of the extent 
of the Stern layer. As shown in Fig. [TJ the cumulative 
counterion charge at this point, q , is significantly smaller 
than the bare charge of the macroion (e.g. qo w —200 e for 
E = 0.01 V/A). Not until r is of the order of 3-4 R M , the 
counterion charge q(r) is equal to the macroion's charge, 
thus completely neutralizing the latter. For high values 
of E, e.g. E — 0.03 V/A, the counterions are stripped off 
the surface of the macroion, since the force due to the 
electric field dominates over the Coulomb attraction be- 
tween macroion and counterions. This leads to a lower 
value of go for high values of E and a less efficient shield- 
ing of the macroion compared to the case of low values of 
E. This in turn increases the velocity Vm of the macroion 
and thus explains the low values of V c t/VM for large val- 
ues of E at distances far from the colloidal surface. 

One might expect, that the region between the Stern 
layer and the point where q ~ Q, is strongly affected 
by hydrodynamic flow features. That this is indeed the 
case can be infered from a comparison to LD simulations 
where the coupling to the LB fluid is switched off (in the 
following, we refer to simulations with a coupling to LB 
as HD simulations). 

In Fig. [21 V c t(r)/VM from LD simulations is compared 
to the same quantity from HD simulations for the two 
macroion charges Q = 255 e and Q = 801 e (in the latter 
case the system contains N ct = 1001 counterions and 
N co = 200 coions). In the HD case, we see that for 
Q = 801 e, the ratio V c t(r)/VM is very close to one in 
the Stern layer region. However, for Q — 255 e the ratio 
is significantly smaller. The strong Coulomb attraction 
dominates the viscous drag and thermal fluctuations for 
high charges, such that layers of counterions nearest to 
the colloidal surface nearly stick to the surface. But for 
lower macroion charges, one starts to see deviations from 
the assumption that the dynamic Stern layer consists of 
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FIG. 2: V c t(r) /Vm as a function of t/Rm for two different 
values of the bare charge Q, as indicated. For the HD case, the 
data from Fig. [T] for Q = 255 e and E = 0.01 V/A are plotted 
(open circles). The closed circles are the corresponding results 
from LD simulations. Also shown are results for Q = 801 
(open and closed diamonds for HD and LD, respectively). In 
this case, the system contains N ct = 1001 counterions and 
N co = 200 coions. 



immobile counterions 0, [IH . Moreover, the Stern layer is 
extended, and it is not restricted to one layer of microions 
closest to the macroion surface as observed in [18( . This 
is due to the much higher values of the surface charge 
density a used in our study. 

As indicated by Fig. [2 in the LD case the motion of 
the counterions is less correlated to that of the macroion 
than in the HD case. The value of V c t(r)/VM is smaller 
in the Stern layer region than in the corresponding HD 
data. Moreover, the counterions reverse their velocity at 
about 0.15 Rm away from the macroion surface, followed 
by a more rapid decrease of V c t(r) /Vm than in HD. From 
data presented below (in Fig. rjj, we will see that the 
number of counterions carried along in the Stern layer is 
almost the same for LD and HD, which reveals that, in 
the LD case, the counterions are more densely packed in 
the Stern layer. 

In order to investigate the mobility \i of the macroion 
as a function of its surface charge density a = Q/(4ttR m ), 
two different types of simulations were performed. First, 
the charge Q was varied from Q = 121 e to 801 e while 
keeping the radii i?H and Rm fixed at 10 A and 20 A, 
respectively. The number of counter- and coions used 
for a given value of Q are listed in Table HI Also in- 
cluded in this list is the Debye screening parameter 
K = [47tAb (N c t+N co ) / L 3 ] - 1 / 2 which was roughly kept con- 
stant around 0.13 A -1 . Secondly, runs for Q = 255 e and 
Q = 401 e were done in which the charge density a was 
varied by choosing different macroion radii 10 A< R M < 
20 A in steps of 2 A or 2.5 A (note that the radius i?n 
remains fixed at 10 A in these runs). The results for the 
different runs are shown in Fig. [3] where /x is plotted as a 
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FIG. 3: The macroion mobility [i as a function of the surface 
charge density a for LD and HD simulations as indicated. The 
charge density a is either varied by changing the radius Kyi of 
the macroion from 10 Ato 20 A keeping Q fixed at Q = 255 e 
or Q = 401 e, or by changing Q from Q = 121 e to Q = 801 e 
keeping Rm fixed at 20 A. The number of counterions and 
coions used for each value of Q is mentioned in text. For all 
data, the "hydrodynamic radius" is chosen to be constant at 
Ru = 2.5 a = 10 A. 
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FIG. 4: Ratio qo/Q as a function of a. As indicated, in the 
different data sets either Q or 7?m are fixed while a is varied 
(see also Fig. [3| . 



function of a for both HD and LD runs. In both cases, the 
mobility fi decreases with increasing a. This behavior is 
in qualitative agreement with theoretical calculations for 

high values of the surface charge, i.e. a > 0.01 e/A (see 
Ref. [HI and references therein). Note that, at small a 
the opposite behavior is observed, i.e. an increase of fi 
with a [H. 

For monovalent microions, Fig. [^demonstrates that at 



constant screening parameter k the mobility // is con- 
trolled by the surface charge density of the macroion. 
Discrepancies that can be seen in the plot might stem 
from the slight variation of k in the different runs (see 
Table [TJ. The LD runs that are shown in the figure ex- 
hibit a similar qualitative behavior. However, compared 
to HD, the LD data is shifted towards lower values of /i. 
This is because the hydrodynamic flow field is switched 
off in the LD simulations. The coupling of the LB fluid 
to the motion of the macroion leads to a backflow ef- 
fect in the LB fluid which enhances the mobility of the 
macroion. This backflow effect yields also stronger cor- 
relations between the motion of the counterions and that 
of the macroion. Thus, as shown in Fig. [51 the reduced 
counterion velocity V c t(r)/VM decays slower in the HD 
case than in the LD case. 

The correlations between counterions and macroions 
are further considered in Fig. [3J Here, qo/Q is plotted 
vs. cr. As before, we vary a by changing Rm keeping Q 
fixed or by changing Q with Rm fixed at 20 A. As we 
see, the different data sets fall roughly onto one master 
curve (note that this is even true for the LD data). The 
functional behavior of qo/Q reflects the one for fi. At 
small values of cr, the ratio qo/Q is "rapidly" increasing 
(associated with a rapid decrease of ji) and it seems to 
saturate at high a (as fj, does). Thus, the electrophoretic 
mobility shows a saturation when qo/Q is approaching 
unity. In the latter case, the electric field sees a "particle" 
with an effective charge Q c ff = Q — qo- 

The regime of low macroion charges Q has been studied 
recently by Lobaskin et al. [l4| using a similar LB/MD 
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FIG. 5: Plot of the measure of distortion d of the Stern layer, 
normalised by the Coulomb radius of colloid Rm as a function 
of the inverse surface charge density 1/a (for the definition of 
'd' see text). The inset shows q^d/QRyi which is a measure 
of the charge separation due to the distortion of the EDL. 
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technique. These authors considered a primitive model 
in the salt free regime (i.e. without coions). Their re- 
sults agree with experimental data [l5| for low Q, where 
an increase of \i with Q (a) is observed. In this regime, 
/1 is controlled by the bare charge Q as there is no sig- 
nificant dynamic Stern layer to create a shielding for Q. 
Hence the peak in [i observed in the experimental data of 
Martin-Molina et al. [l3[ denotes the point of crossover 
from this regime to a regime where [i is controlled by 
Qo/Q- For high values of ct, /i decreases with a due to 
greater shielding for higher Q, whereas, for low values of 
a (Q) the mobility increases with a. 

O'Brien and White [7j have proposed that the decrease 
in macroion mobility can be attributed to the increas- 
ing distortion of the EDL at high values of a. To see 
whether the distortion of the charge cloud is indeed a 
relevant effect also at very high values of a (considered 
in this work), we quantify the distortion of the EDL as 
follows: We compute the distance d of the center of mass 
(CM) of the EDL from the center of the macroion. Here, 
those counterions form the EDL which are within the dy- 
namic Stern layer and move along the macroion. Figure[5] 
shows that the normalised distortion dj Rm and the nor- 
malised quantity qod/QRn decrease with increasing a. 
The quantity q^d/QRui is a measure of the charge sep- 
aration between the center of the macroion and center 
of the counterion cloud. Thus, the amount of distortion 
becomes less pronounced with increasing a. This rules 
out the mechanism proposed by O'Brien and White in 
the case of high values of a. 

We have already infered from Fig. [5] that hydrody- 
namic backflow due to the LB fluid enhances the correla- 
tions between the counterion and the macroion motion as 
well as the absolute value of p.. In Fig. [SJ the mobilities 
obtained from LD and HD runs (denoted by /^hd and 
/ild, respectively) are directly compared to each other 



by plotting the ratio /xhd/a'ld as a function of a. We see 
that /^hd/Vld increases almost linearly with a. 

Microscopically, we understand this by noting that the 
higher macroion charge density leads to a higher num- 
ber density of counterions per unit volume in the mov- 
ing EDL of the colloid. The counterion coupling to 
the LB fluid results in an increased velocity of the LB 
fluid around the macroion at higher densities. Thus the 
macroion feels less frictional drag and, compared to the 
LD case, a higher value of is obtained. 

In Fig. [7J it is further explored how the counterion 
coupling with the LB fluid affects the colloidal mobility. 
We keep the friction coefficient for the macroion coupling 
to the LB fluid, £0, fixed and vary that for the microion 
coupling, £06 (for the definition see Sec. II). Whereas, for 
the HD simulations, fi decreases significantly with £ofei in 
the LD case, /j, exhibits only a weak dependence on £ofc- 
Thereby, the HD result tends to approach the one from 
LD for high values of £ob- In the HD case, small values 
of £ot> mean a weaker coupling of the microion motion 
to that of the LB fluid, leading to an increase of the 
macroion mobility. The detailed understanding of this 
finding requires further investigation. 

Note that, in the LD simulation for Q — 121 e, a 
dynamic Stern layer cannot be formed, because the 
Coulomb attraction is too weak in this case. Hence, a 
high value of mobility results which is outside the range 
used in Fig. [7] On the other hand, in the HD simulation 
for Q = 121 e a well-defined dynamic Stern layer is seen 
and, as shown by Fig.[7l the behavior of /i is qualitatively 
similar to that at Q = 255 e. This demonstrates the im- 
portance of the hydrodynamic flow field for the behavior 
of the electrophoretic mobility of weakly charged colloids. 
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radius used. 



ulations is compared to that from LD simulations. 

To this end, two systems of a macroion of charge 
Q = 255 e mixed with divalent microions are considered, 
one with N ct — 526 divalent counterions and 7V co = 400 
divalent coions, and the other one with N ct = 176 di- 
valent counterions and N co = 50 divalent coions. In 
each case, 3 monovalent counterions are added to main- 
tain charge neutrality. The charge density a is varied 
by changing the radius of the macroion Rm from 10 to 
20 A, thereby keeping the radius i?n fixed at 10 A. In 
Fig. [51 the mobility fi as function of a is plotted for the 
different systems. Different from the case of monovalent 
microions, fj, changes its sign in all cases, thus indicating 
that at high values of a the macroion is moving opposite 
to the direction of the electric field. Overall, at a given 
value of <7, the mobility of the macroion in the presence of 
divalent microions is smaller than with monovalent ones. 

Since the screening due to divalent microions is much 
more effective, the counterion charge density is higher 
near the colloidal surface than in electrolytes with mono- 
valent microions. This can be seen in Fig. [5J where the 
system with iV co = 400 coions is considered. The cu- 
mulative charge q(r) around the macroion is displayed 
for different values of Rm corresponding to different val- 
ues of a at the fixed macroion charge Q = 255 e. We 
see that for Rm < 16 A, q(r) develops a local minimum 
near the colloidal surface. This indicates charge inver- 
sion (i.e. \q(r)\ > 255 e) and it explains the appearance 
of negative values of fi. 

For comparison, we have plotted LD data in Fig. [5] 
for the largest and the smallest radii used, i.e., for Rm = 
20 A and Rm = 10 A, respectively. We see that the charge 
profile for LD and HD is nearly identical, which shows 
that the dynamic Stern layer is determined by the static 
Coulomb forces. For higher salt concentration, charge 
inversion is more pronounced, as seen by lower values of 
/j, at high a. Both LD and HD show the phenomenon of 
charge inversion at high values of a. In particular in the 
HD case, charge inversion is more pronounced for higher 
salt concentration. As in the studies with monovalent 
counter- and coions, for divalent microions | /xld |<| 
Mhd I holds as well. 



IV. CONCLUSIONS 



B. Divalent salt ions 

Charge inversion in the presence of multivalent coun- 
terions has been observed in static experiments and the- 
oretical studies [HI, HH, [36| . Experimentally charge inver- 
sion is detected using mainly electrophoresis. It is well 
established that charge inversion is due to counterion cor- 
relations and beyond the scope of mean-field theories. 
Our aim in this study is twofold: the onset charge den- 
sity cr, at which charge inversion occurs, is determined 
and secondly, the extent of the Stern layer from HD sim- 



In this work, a hybrid MD/LB scheme was used to in- 
vestigate the electrophoretic mobility of a macroion in an 
electrolyte solution. We have considered highly-charged 
macroions (a = 38 /zC/cm 2 -500 /xC/cm 2 ) for which the- 
ories based on the linearized electrokinetic equation are 
not applicable. For high values of the macroion's surface 
char ge g , experiments on systems with multivalent salt 
ions [13| have shown that the electrophoretic mobility fi 
can decrease with increasing a for 2 : 1 salts whereas fi 
for 2 : 2 salts, it shows a plateau. In our simulations, we 
observe decrease of fj, with a for both 1 : 1 and 2 : 2 salt. 
Moreover, in contrast to previous theoretical studies 0], 
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we observe a lowering of /j, with a even for kRm < 3 for 
monovalent salt ions. 

The authors of [HI , attribute the lowering of \x in the 
2 : 1 case to the presence of coions in the EDL. As we have 
demonstrated in this work, a decreasing \i as a function of 
a requires the formation of a pronounced dynamic Stern 
layer of charge q$, consisting of counterions that move 
along the same direction as the macroion in the presence 
of an electric field. The lowering of fj, is directly related 
to the increase of the quantity qo/Q which we call the 
screening charge fraction. The absolute value of screen- 
ing charge fraction \qo/Q\ of the colloid-counterion cloud 
complex increases with increasing cr, reflecting a lower- 
ing of fj,. Note that the number of coions in the dynamic 
Stern layer is negligible. 

Furthermore, we observe charge inversion for divalent 
counter- and coions, but only when we reach sufficiently 
high values of a. At very high values of cr, /i for colloidal 
systems with monovalent salt ions tends to saturate, be- 
cause there is no space to add more counterions into the 
dynamic Stern layer. However, this is different when 
multivalent salt ions and counterions are used. Then, 
l?o| IQ > 1 holds and \i becomes negative. The distribu- 
tion of charge in the EDL for this case as a function of dis- 
tance r from the center of the colloid is non-monotic and 
completely different from what is expected from standard 
electrokinetic theories [5j. Though we have not carried 
out studies for the 2:1 salt case, it is possible that the dis- 
crepancy in the experimental results between the 2:1 and 
the 2:2 salt case arises from entropic/osmotic forces. The 
larger number of coions in the 2:1 salt (compared to the 
systems with a 2:2 salt) may lead to larger osmotic forces 
on the counterions which could result in larger values of 
qo for the 2:1 case. However, this issue can be resolved 
only after further investigations. 

In the work of O'Brien and White Q, the lowering of 



/i with the C potential has been attributed to the dis- 
tortion of the EDL. Though the explicit measurement 
of the retarding force arising from the distortion of the 
EDL is outside the scope of this work, we have shown 
that the charge distortion within the dynamic Stern layer 
decreases with increasing cr. This is indicated by the low- 
ering of the quantities d/Ryi and qod/QRM with increas- 
ing a. This observation combined with a host of other 
evidences presented in this paper, indicate that accumu- 
lation of charges in the Stern layer is the dominant mech- 
anism for the lowering of fi as a function of cr, providing 
that high values of cr are considered. 

By comparing the LB/MD (or HD) results to those 
from LD simulations, we were able to elucidate the role 
of hydrodynamic interactions. The structure of the Stern 
layer is more compact in the Langevin simulations. Hy- 
drodynamic interactions enhance the electrophoretic mo- 
bility and so they have the opposite effect to electrostatic 
screening. This is due to a backflow effect in the fluid 
which pushes the counterions to move along the same di- 
rection as the macroion. Hydrodynamic forces have pro- 
nounced influence in the formation of the dynamic Stern 
layer, especially for low charges. An interesting finding 
is that the mobility as obtained from HD relative to that 
from LD exhibits an almost linear increase with cr. Fur- 
ther work on the role of hydrodynamic interactions for 
the electrophoresis of charged colloids is in progress. 
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